#include<iostream>
#include<stdio.h>
#include<stdlib.h>
#include <iomanip>
#include <sys/time.h>
#include"cuda_runtime.h"
#include"device_launch_parameters.h"
using namespace std;
const int BLOCK_SIZE=1024;
float *A;
int n=500;
void normal(float* a)
{
   
    int i, j, k;
    for (k = 0; k < n; k++)
    {
        for (j = k + 1; j < n; j++)
        {
            a[k*n+j] = a[k*n+j] / a[k*n+k];
        }
        a[k*n+k] = 1.0;
        for (i = k + 1; i < n; i++)
        {
            for (j = k + 1; j < n; j++)
            {
                a[i*n+j] = a[i*n+j] - a[i*n+k] * a[k*n+j];
            }
            a[i*n+k] = 0;
        }
    }
}
void m_rest(float* A, int n)
{
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < i; j++)
            A[i*n+j] = 0;
        A[i*n+i] = 1.0;
        for (int j = i + 1; j < n; j++)
        {
            float x = rand();
            while (x == 0) x = rand();
            A[i*n+j] = x;
        }
    }
    for (int k = 0; k < n; k++)
        for (int i = k + 1; i < n; i++)
            for (int j = 0; j < i; j++)
                A[i*n+j] += A[k*n+j];
}

__global__ void division_kernel(float* data, int k, int n)
{
    int tid = blockDim.x * blockIdx.x + threadIdx.x;
    int element = data[k*n+k];
    int temp = data[k*n+tid];
    data[k*n+tid] = (float)temp/element;
    return;  
}
__global__ void eliminate_kernel(float* data, int k, int n){
     int tx = blockDim.x * blockIdx.x + threadIdx.x;
     if(tx==0)  data[k*n+k]=1.0;//对角线元素设为 1 
    int row = k+1+blockIdx.x;//每个块负责一行
     while(row<n){
     int tid = threadIdx.x;
     while(k+1+tid < n){
     int col = k+1+tid;
     float temp_1 = data[(row*n) + col];
     float temp_2 = data[(row*n)+k];
     float temp_3 = data[k*n+col];
    data[(row*n) + col] = temp_1 - temp_2*temp_3;
 tid = tid + blockDim.x;
 }
 __syncthreads();//块内同步
 if (threadIdx.x == 0){
 data[row * n + k] = 0;
 }
 row += gridDim.x;
 }
 return;
}
int main()
{
    struct timeval tstart, tend;
    double timeUsed0=0.0;
    A = new float[n*n];
    float*B = new float[n*n];
    float* gpudata;
    m_rest(A, n);
    m_rest(B, n);
    gettimeofday(&tstart, NULL);
    normal(B);
    gettimeofday(&tend, NULL);
    timeUsed0 = 1000000 * (tend.tv_sec - tstart.tv_sec) + tend.tv_usec - tstart.tv_usec;
    dim3 dimBlock(BLOCK_SIZE,1);
    dim3 dimGrid(1, 1);
    cudaEvent_t start, stop;
    float elapsedTime = 0.0;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    cudaError_t ret;
    ret=cudaMalloc(&gpudata,n*n*sizeof(float));
    ret=cudaMemcpy(gpudata,A,n*n,cudaMemcpyHostToDevice);
    cudaEventRecord(start, 0);
    for(int k=0;k<n;k++){
    division_kernel<<<dimGrid,dimBlock>>>(gpudata,k,n);//负责除法任务的核函数
    cudaDeviceSynchronize();//CPU 与 GPU 之间的同步函数
    ret = cudaGetLastError();
    if(ret!=cudaSuccess){
    printf("division_kernel failed, %s\n",cudaGetErrorString(ret));
    } 
    eliminate_kernel<<<dimGrid,dimBlock>>>(gpudata,k,n);//负责消去任务的核函数
    cudaDeviceSynchronize();
    ret = cudaGetLastError();
    if(ret!=cudaSuccess){
    printf("eliminate_kernel failed, %s\n",cudaGetErrorString(ret));
    }
    }
    cudaEventRecord(stop, 0);
    ret=cudaMemcpy(A, gpudata, n*n, cudaMemcpyDeviceToHost);
    cudaEventSynchronize(stop);
    cudaEventElapsedTime(&elapsedTime, start, stop);
   printf("CPU_LU:%f ms\n",  timeUsed0/1000);
    printf("GPU_LU:%f ms\n", elapsedTime);
    cudaEventDestroy(start);
    cudaEventDestroy(stop);
}